#### Import necessary modules

import os

import pandas as pd

from pathlib import Path

import copy

import jwrangle

import jvis

import jinspect

import jtest

import jweb

from imp import reload

import upsetplot as upset

from Bio import SeqIO

import altair as alt

​

import numpy as np

import preprocess

import seaborn as sns

​

#### define working directories

cwd  = Path(os.getcwd())

base_path = Path(os.path.join(*cwd.parts[:cwd.parts.index('experiments')]))

​

#### MaxQuant proteinGroups & evidence files

MQ_folder = jwrangle.importMixedFiles(cwd / 'MaxQuant')

MQ_folder.keys()

​

pGroups = MQ_folder['proteinGroups.txt']

evidence = MQ_folder['evidence.txt']

​

#### Inspect MQ setup

MQ_folder['parameters.txt'].head(9)

#### Inspect column names

colnames = list(pGroups.columns.values)

​

#### Derive experiment names as a list

#### I gave the original samples some irritating names that I'd like to correct

experiment_names = []

for i in colnames:

    if 'Intensity ' in i:

        experiment_names.append(i.replace('Intensity ', ''))

​

new_expt_names = []

for i in experiment_names:

    i = i.replace('pH8S', 'pH7-S')

    i = i.replace('ne', 'ne-')

    i = i.replace('NC', 'nCL')

    i = i.replace('PC', '254')

    i = i.replace('heat', '-heat')

    i = i.replace('obd', '-obd')

    new_expt_names.append(i)

#### Create a list of associated conditions

#### These samples were mislabeled 'NC' and 'PC' (the latter is often short for PAR-CL, but these samples are 254nm crosslinked)

condition =  ['b3e3_nCL']*4 + ['b3e3_254']*4 + \

             ['b4e3_nCL']*4 + ['b4e3_254']*4 + \

             ['b3e4_nCL']*4 + ['b3e4_254']*4 + \

             ['b4e4_nCL']*4 + ['b4e4_254']*4

​

​

#### Create a list of associated replicate identifiers

replicate = []

for i in new_expt_names:

    replicate.append(i.split('-')[-1])

​

#### Create a more reader friendly list of sample names

samples =  []

for pos, label in enumerate(condition):  

    samples.append(label+new_expt_names[pos][-2:])

#### The MQ_groups dictionary assigns each condition to its MQ group parameter

MQ_groups = {'b3e3':['b3e3_nCL','b3e3_254'], 

             'b4e3':['b4e3_nCL','b4e3_254'],

             'b3e4':['b3e4_nCL','b3e4_254'],

             'b4e4':['b4e4_nCL','b4e4_254']}

​

#### We then use the condition list to create an MQ_groups list

Grp_parameters = []

for label in condition:

    for key, value in MQ_groups.items():

        if label in value:

            Grp_parameters.append(key)

​

exp_df = pd.DataFrame(

    {'experiment': experiment_names,

     'condition': condition,

     'replicate': replicate,

     'sample':samples,

     'measure':['Intensity']*len(samples),                  # adding this column allows our metadata file to be compatible with Proteus

     'MQgroups':Grp_parameters

    })

​

exp_df

MQ_expt296 = jwrangle.MQ_writeMetadata(pGroups, evidence, exp_df, 'e296r', cwd)

#### If not already loaded, read in the metadata-adjusted files

metadata = pd.read_csv(cwd / 'metadata' / 'e296r_metadata.csv', index_col = 0)

pGroups = pd.read_csv(cwd / 'MaxQuant' / 'e296r_proteinGroups_metalabeled.txt', delimiter = '\t')

evidence = pd.read_csv(cwd / 'MaxQuant' / 'e296r_evidence_metalabeled.txt', delimiter = '\t')

#### If not already loaded, read in the remapped proteinGroups file

pGroups_remap = jwrangle.importMixedFiles(cwd / 'downloads' / 'pGroups_remap', dropSuffix = 'yes')

pGroups_remap.keys()

#### jwrangle.importMixedFiles() returns a dictionary where keys = files. We want the 'id_map' table created by jweb.mapAnyID_gPro().

#### We'll rename the Query column and drop duplicates so the table can be merged with our proteinGroups table.

id_map = pGroups_remap['id_map'].rename(columns={'Query': 'Majority protein IDs'}).drop_duplicates()

​

#### Now use merge to add these new columns to our proteinGroups table

pGroups_map = pd.merge(pGroups, id_map, on='Majority protein IDs', how='left')

​

#### Check the tables are merged by viewing column elements from each.

pGroups_map[id_map.columns.tolist() + ['Peptide IDs']].head(2)

#### Extract contaminants

contaminants = jinspect.MQ_getContaminants(pGroups_map, metadata)

​

#### Sort the metadata into a more intuitive order

metadata_sort = metadata.sort_values(by = ['MQgroups','condition', 'replicate'], ascending=[True, False, True])

​

#### Visually inspect contaminants  

jvis.MQ_getContaminants_sbplot(contaminants, metadata_sort, width = 3, length = 2, layout = 'single')

#### Extract the missed cleavage data into a long form table for plotting

MissedCleavages = jinspect.MQ_getMissedCleavages(evidence, metadata, drop_contaminants = True)

MissedCleavages.sort_values(by=['expt','sample'], inplace = True)

MissedCleavages

#### Map the location of the custom FASTA elements

os.listdir(base_path / 'my_resources' / 'FASTA')  

​

#### Create a list of the non-human proteins that were added to the custom FASTA genome search. 

new_cont = []

with open(base_path / 'my_resources' / 'FASTA' / "custom_proteome_elements.fasta", "r") as handle:

    for record in SeqIO.parse(handle, "fasta"):

        new_cont.append(record.id.split('|')[1])

​

#### Remove all unwanted contaminants and IDs from the proteinGroups table      

pGroup_clean = jwrangle.MQ_getThreePassFilter(pGroups_map, custom_exclusion = new_cont)

​

#### Inspect the cleaned dataframe

pGroup_clean[['ENTREZGENE_gPro primary'] + [i for i in pGroup_clean.columns if 'iBAQ' in i]].head(2)

#### Drop duplicates and apply LFQ filter

filter_dict = jinspect.MQ_dropDuplicateIDs(pGroup_clean, metadata, prefix = 'Peptides', ID = 'ENTREZGENE_gPro primary', pool = 'measure', drop_ID = 'None', 

                                            keep_PoolCalcs = False, applyLFQ_filter = ['Intensity', 'iBAQ'])

​

#### The df_keep value contains our targets, df_droprows conatins the discarded duplicates. Assign the df_keep value to a new variable and inspect.

pGroup_filtered = filter_dict['df_keep']

pGroup_filtered.head(2)

#### Confirm that clustering matches expectations

jvis.MQ_showDendrogramQC_mplplot(pGroup_filtered, 'Intensity', metadata, 'QC clustering: ', grid = 'YES', fsize = (8, 8))

#### Log2 transform available intensity values.

pGroup_log2 = jwrangle.Log2_ByPrefix(pGroup_filtered, 'LFQ intensity')

pGroup_log2 = jwrangle.Log2_ByPrefix(pGroup_log2, 'iBAQ')

pGroup_log2 = jwrangle.Log2_ByPrefix(pGroup_log2, 'Intensity')

pGroup_log2.replace(0,np.nan, inplace=True)

​

#### Create a long form dataset for each desired grouping

pool_SampleIntensity = jwrangle.MQ_poolMulti(pGroup_log2, metadata, melt_list = ['Intensity', 'LFQ intensity'], group = 'condition')

pool_SampleIntensity.keys()

#### Inspect the Intensity shifts generated by the LFQ normalisation algorithm. In MaxQuant Raw Intensity and iBAQ values are 

#### not subjected to the LFQ normalisation calculations so this is a good way to spot any gross violations

sns.set_style('whitegrid')

jvis.ViolinCompare_sbplot(pool_SampleIntensity['b3e3_nCL'], title = 'b3e3_nCL: Normalisation Effects', ylabel = 'Log2(Intensity)', palette = ['#ff6666', '#99ccff'])

#### Inspect the Intensity shifts generated by the LFQ normalisation algorithm. In MaxQuant Raw Intensity and iBAQ values are 

#### not subjected to the LFQ normalisation calculations so this is a good way to spot any gross violations

sns.set_style('whitegrid')

jvis.ViolinCompare_sbplot(pool_SampleIntensity['b3e3_254'], title = 'b3e3_254: Normalisation Effects', ylabel = 'Log2(Intensity)', palette = ['#ff6666', '#99ccff'])

#### Inspect the Intensity shifts generated by the LFQ normalisation algorithm. In MaxQuant Raw Intensity and iBAQ values are 

#### not subjected to the LFQ normalisation calculations so this is a good way to spot any gross violations

sns.set_style('whitegrid')

jvis.ViolinCompare_sbplot(pool_SampleIntensity['b3e4_nCL'], title = 'b3e4_nCL: Normalisation Effects', ylabel = 'Log2(Intensity)', palette = cpal['Set3_qual'])

#### Inspect the Intensity shifts generated by the LFQ normalisation algorithm. In MaxQuant Raw Intensity and iBAQ values are 

#### not subjected to the LFQ normalisation calculations so this is a good way to spot any gross violations

sns.set_style('whitegrid')

jvis.ViolinCompare_sbplot(pool_SampleIntensity['b3e4_254'], title = 'b3e4_254: Normalisation Effects', ylabel = 'Log2(Intensity)', palette = cpal['Set3_qual'])

#### Inspect the Intensity shifts generated by the LFQ normalisation algorithm. In MaxQuant Raw Intensity and iBAQ values are 

#### not subjected to the LFQ normalisation calculations so this is a good way to spot any gross violations

sns.set_style('whitegrid')

jvis.ViolinCompare_sbplot(pool_SampleIntensity['b4e3_nCL'], title = 'b4e3_nCL: Normalisation Effects', ylabel = 'Log2(Intensity)', palette = cpal['Set3_qual'][4:])

#### Inspect the Intensity shifts generated by the LFQ normalisation algorithm. In MaxQuant Raw Intensity and iBAQ values are 

#### not subjected to the LFQ normalisation calculations so this is a good way to spot any gross violations

sns.set_style('whitegrid')

jvis.ViolinCompare_sbplot(pool_SampleIntensity['b4e3_254'], title = 'b4e3_254: Normalisation Effects', ylabel = 'Log2(Intensity)', palette = cpal['Set3_qual'][4:])

#### Inspect the Intensity shifts generated by the LFQ normalisation algorithm. In MaxQuant Raw Intensity and iBAQ values are 

#### not subjected to the LFQ normalisation calculations so this is a good way to spot any gross violations

sns.set_style('whitegrid')

jvis.ViolinCompare_sbplot(pool_SampleIntensity['b4e4_nCL'], title = 'b4e4_nCL: Normalisation Effects', ylabel = 'Log2(Intensity)', palette = cpal['Set1_qual'][2:])

#### Inspect the Intensity shifts generated by the LFQ normalisation algorithm. In MaxQuant Raw Intensity and iBAQ values are 

#### not subjected to the LFQ normalisation calculations so this is a good way to spot any gross violations

sns.set_style('whitegrid')

jvis.ViolinCompare_sbplot(pool_SampleIntensity['b4e4_nCL'], title = 'b4e4_nCL: Normalisation Effects', ylabel = 'Log2(Intensity)', palette = cpal['Set1_qual'][2:])

#### Pool data into a single long form dataset

pooled_dfDropGroupOne = jwrangle.MQ_poolDataByCondition(pGroup_log2, metadata_sort, prefix_list = ['Intensity', 'Sequence coverage'])

​

#### Compare Intensity distribution using a box and whisker plot

sns.set_style('whitegrid')

jvis.BoxPlotByColumn_sbplot(pooled_dfDropGroupOne, 'Intensity: ', 'Intensity')

#### Compare Sequence coverage using a box and whisker plot

sns.set_style('whitegrid')

jvis.BoxPlotByColumn_sbplot(pooled_dfDropGroupOne, 'Sequence coverage: ', 'Sequence coverage %')

#### Extract the total peptides observed per sample

metaStats = jinspect.MQ_getSumBySample(pGroup_log2, metadata_sort, freqList = ['Peptides'], measure = False)

​

#### Plot the sum peptides

sns.set_style('whitegrid')

jvis.BarPlotByGroup_sbplot(metaStats, x_col = 'condition', y_col = 'Peptides', title = 'Sum Peptides vs pH enrichment', pal = set2_paired,

                          errorbars = 'SEM')

#### Count the number of unique 

metaStats = jinspect.MQ_getFrequencyBySample(pGroup_log2, metaStats, freqList = ['Intensity', 'iBAQ', 'LFQ intensity'], measure = False)

metaStats.iloc[:,1:]

#### Plot the counts

sns.set_style('whitegrid')

jvis.BarPlotByGroup_sbplot(metaStats, x_col = 'condition', y_col = 'Intensity', title = '# Genes Detected By Group', pal = set2_paired, ylabel = 'Unique Genes', errorbars = 'SEM')

#### Extract the intensity values as a dictionary where keys = groups

Intensity_Dict = jwrangle.MQ_getSliceByPrefix(pGroup_log2, metadata, 'Intensity', group = 'condition', add_col = None)

Intensity_Dict.keys()

#### Use the metadata and proteinGroups tables to count how many times a gene is identified in its group (/6). 

#### Here I demonstrate how we can count for all instances of Intensity, iBAQ and LFQ Intensity

pGroup_Freq = jinspect.MQ_getFrequencyByGroup(pGroup_log2, metadata, 'iBAQ', group = 'condition')

pGroup_Freq = jinspect.MQ_getFrequencyByGroup(pGroup_Freq, metadata, 'LFQ intensity', group = 'condition')

pGroup_Freq = jinspect.MQ_getFrequencyByGroup(pGroup_Freq, metadata, 'Intensity', group = 'condition')

pGroup_Freq[['ENTREZGENE_gPro primary'] + [i for i in pGroup_Freq.columns if 'Freq' in i]].head(2)

#### Now we can classify RBP; the console will report if all proteins/genes have been properly classified or not.

RBP_Dict_b3e3 = jtest.MQ_applyClassifyRBP(pGroup_Freq, 'LFQ intensity', metadata, 'b3e3', 'b3e3_nCL', 'b3e3_254', 3, 

                                         add_cols = ['ENTREZGENE_gPro all', 'ENTREZGENE_gPro primary', 'ENTREZGENE_gPro name'])

#### Now we can classify RBP; the console will report if all proteins/genes have been properly classified or not.

RBP_Dict_b4e3 = jtest.MQ_applyClassifyRBP(pGroup_Freq, 'LFQ intensity', metadata, 'b4e3', 'b4e3_nCL', 'b4e3_254', 3, 

                                         add_cols = ['ENTREZGENE_gPro all', 'ENTREZGENE_gPro primary', 'ENTREZGENE_gPro name'])

#### Now we can classify RBP; the console will report if all proteins/genes have been properly classified or not.

RBP_Dict_b4e4 = jtest.MQ_applyClassifyRBP(pGroup_Freq, 'LFQ intensity', metadata, 'b4e4', 'b4e4_nCL', 'b4e4_254', 3, 

                                         add_cols = ['ENTREZGENE_gPro all', 'ENTREZGENE_gPro primary', 'ENTREZGENE_gPro name'])

#### The results of this classifcation can be found in a dictionary of dataframes for each output

#### See help() for an explanation of each dataset

RBP_Dict_b3e3.keys()

#### Because each subclass of the class II RBP are identified based on different statistical assumptions we can more closely inspect if those assumptions trend differently among the different treatments.

plt.figure('rbp class')

sns.set_style('whitegrid')

sns.set(font_scale=1.1)

ax = sns.countplot(x='MQgroup', hue = 'RBP subClass', data=RBP_Class[RBP_Class['RBP subClass']!=''].sort_values(by=['MQgroup','RBP subClass']), palette = cpal['Set3_qual'], edgecolor = 'black')

ax.set_ylabel('Unique Gene Count')

ax.set_title('Unique Genes detected per RBP Class')

​

RBP_Class.head(1)

from matplotlib import pyplot as plt

import numpy as np

from matplotlib_venn import venn3, venn3_circles

​

plt.figure(figsize=(6,6))

​

v = venn3([set(RBP_Dict_b3e4['NC']['ENTREZGENE_gPro primary']), 

           set(RBP_Dict_b3e3['I']['ENTREZGENE_gPro primary']),

           set(RBP_Dict_b4e3['I']['ENTREZGENE_gPro primary'])], 

           set_labels = ('b3e4', 'b3e3', 'b4e3'),

           alpha = 0.5)

​

c = venn3_circles([set(RBP_Dict_b3e4['NC']['ENTREZGENE_gPro primary']), 

                   set(RBP_Dict_b3e3['I']['ENTREZGENE_gPro primary']), 

                   set(RBP_Dict_b4e3['I']['ENTREZGENE_gPro primary'])], linestyle='dashed')

​

c[0].set_lw(1.5)

c[0].set_ls('dotted')

c[1].set_lw(1.5)

c[1].set_ls('dotted')

c[2].set_lw(1.5)

c[2].set_ls('dotted')

​

for text in v.set_labels:

    text.set_fontsize(18)

for text in v.subset_labels:

    text.set_fontsize(14)

#### I like to keep a table of interesting GO terms in a local csv file, let's find it

MyResources = jwrangle.importMixedFiles(base_path / 'my_resources')

MyResources.keys()

MyResources['GO_TermsOfInterest.csv'].head(9)

#### Let's get basic statistics on the broader RNA-binding category "GO:0003723"  

MyGO_Stats_Dict = jweb.fetchQuickGO_stats(['GO:0003723'], QG_geneProductType = 'protein', QG_taxonId = '9606', QG_geneProductSubset = ['Swiss-Prot', 'TrEMBL'])

MyGO_Stats_Dict.keys()

#### To view basic information, like the number of annotations associated with the GO term

RNA_Binding_stats = jweb.getQuickGO_stats(MyGO_Stats_Dict)

RNA_Binding_stats

#### Download and create a local copy of all the records associated with our GO query.  

QuickGo_dict = jweb.fetchQuickGO(['GO:0003723'], QG_geneProductType = 'protein', QG_taxonId = '9606', QG_geneProductSubset = ['Swiss-Prot', 'TrEMBL'], gConvertOrganism='hsapiens', writetopath= True)

QuickGo_dict.keys()

### We'll assess our RBP frequency by considering both SwissProt and TrEmbl records, thus let's concatenate the dictionaries for GO:0003723

QuickGo_dict_concat = jwrangle.concatGO_DataFrameDict(QuickGo_dict)

QuickGo_dict_concat.keys()

#### With our GO records in hand we next investigate whether our identified proteins are members of GO:OOO3723

#### This can be done with any DataFrame, but here we'll use the last version of our modified proteinGroups table 

pGroup_GO = jwrangle.AnnotateDataFrameCtrls(pGroup_Freq, QuickGo_dict_concat, search_match = 'ENTREZGENE_gPro primary', dict_match = 'ENTREZGENE_gPro primary', none_col = 'GO_None')